////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#include "headers.hh"
#include "BaseUtils.h"
#include "dEdxGausN.h"

//read from tree "mn2d" in outmn2d_post_iBeta100_save.root (iBeta100 is the final result in guided convergence)
//see 
//void MN2D::Save(const TString tag)
//void dEdxGausN::SavePar(TTreeSRedirector * stream, const TString tn,  Double_t *totalyield,  Double_t *yields, Double_t *pars,  Double_t *lowerrs,  Double_t *higherrs, const Int_t order[], Double_t *outpar, Double_t *outlowerr, Double_t *outhigherr, Int_t *orderinv)
//void dEdxGausN::ToGausPar(const Double_t *prepar, Double_t *gauspar, const Double_t *prelowerr, const Double_t *prehigherr, Double_t *gauslowerr, Double_t *gaushigherr)
//void MN2D::SavePhysics(TTreeSRedirector * stream, const TString tn, const Int_t ibin)

/*
totn: real ncount in slice ix
rawp$(particle)$(internal-par): internal parameter (fraction, relative position, resolution) for each particle

y$(particle): yield of each particle = totn * fraction
p$(particle)$(gaus-par): gaussian parameter (fraction, absolute position, resolution) for each partile 

aa$(particle): fraction, identical to p$(particle)0
cova$(particle): convariance matrix between particle and pion
rr$(particle): ratio to pion

ix: ith bin
xx: log10(p)
 */

void readiBeta(const Int_t ntype, const TString outname)
{
  TList *ll = new TList;

  //totn from horiginalyield
  TGraphAsymmErrors *gtotn=BaseUtils::TreeToGraph("mn2d","1",(Char_t*)"xx",(Char_t*)"totn", (Char_t*)"",(Char_t*)"",(Char_t*)"",(Char_t*)"", kTRUE); 
  gtotn->SetName("TotalYield");
  gtotn->SetTitle("TotalYield vs. momentum");
  ll->Add(gtotn);

  //p0: fraction, p2: resolution, p3 =rr(ratio), all finally as a function of momentum
  const Int_t vidIn[]={0,2,3};
  const TString vnOut[]={"Fraction", "Resolution", "Ratio"};
  const Int_t nvar=sizeof(vnOut)/sizeof(TString);
  TGraphAsymmErrors * gout[nvar][ntype];

  for(Int_t itype=0; itype<ntype; itype++){
    const Char_t * tn = dEdxGausN::ParticleName(itype);
    for(Int_t ivar=0;ivar<nvar; ivar++){
      if(vidIn[ivar]==3){
        gout[ivar][itype]=BaseUtils::TreeToGraph("mn2d","1",(Char_t*)"xx",Form("rr%d", itype),   (Char_t*)"",(Char_t*)"",Form("er%d", itype),Form("er%d", itype), kTRUE); 
      }
      else{
        const Char_t * tag=Form("%d%d", itype, vidIn[ivar]);
        gout[ivar][itype]=BaseUtils::TreeToGraph("mn2d","1",(Char_t*)"xx",Form("p%s", tag),   (Char_t*)"",(Char_t*)"",Form("le%s", tag),Form("he%s", tag), kTRUE); 
      }
      gout[ivar][itype]->SetName(Form("%s%s", tn, vnOut[ivar].Data()));
      gout[ivar][itype]->SetTitle(Form("%s vs. momentum", gout[ivar][itype]->GetName()));
      gout[ivar][itype]->SetMaximum(1.1);
      gout[ivar][itype]->SetMinimum(0);

      ll->Add(gout[ivar][itype]);
    }
  }

  TFile *fout=new TFile(Form("readOutput_%s.root", outname.Data()),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("readiBeta argc!=2\n"); 
    return 1;
  }
  dEdxGausN::Ini();

  //=====

  const TString fin(argv[1]);
  TFile::Open(fin);
  TString outname(fin);
  outname.ReplaceAll("out","");
  outname.ReplaceAll("_post","");
  outname.ReplaceAll("_save","");
  outname.ReplaceAll(".root","");
  readiBeta(dEdxGausN::Ntype(), outname);
  return 0;
}
